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Main  Body 

1  On  the  left  we  show  the  PFBA  (black  circle)  amid  southern  California  topography  2 
and  faults.  The  map  on  the  lower  right  shows  the  PFBA  and  ANZA  network 
stations  (black  triangles).  The  PFBA  (upper  right)  is  a  28  element  6  km  aperture 
array  which  was  deployed  for  3  months  in  early  1991.  Each  element  consists  of  3 
orthogonal  STS-2  sensors  (.0083  to  40  Hz).  The  Salton  Trough  lies  to  the  north 
and  east  of  these  deployments. 

2  Record  section  of  events  recorded  by  the  vertical  component  of  the  ANZA  3 
network  station  PFO  (located  at  the  intersection  of  the  crosshairs  in  the  map).  The 
events  occurred  between  0.71  and  0.78°  from  the  station.  This  record  section 
includes  events  at  depths  ranging  from  0.8  to  5  km.  Ml  ranges  from  1.9  to  3.1. 
Despite  the  similar  propagation  lengths  these  seismic  records  display  little  in 
common. 

3  In  this  figure  we  display  vertical  component  recordings  of  Landers  aftershocks  4 
(above  Ml  3.5)  recorded  by  three  of  the  ANZA  network  stations.  The  Pinon  Flat 
station  shows  that  an  obvious  phase  (delayed  2  s  after  the  onset)  appears  shortly 
before  .7°  and  suggests  that  it  disappears  beyond  0.85°.  The  sudden  appearance 

of  this  phase  suggests  that  it  is  a  bounce  from  the  Moho  that  reaches  the  critical 
angle  by  0.7°.  The  other  two  stations  suggest  that  the  energy  is  not  due  to  a 
simple  bounce  since  there  is  little  evidence  of  a  late  phase  at  these  loeations. 

4  Vertical  component  triggered-mode  recordings  of  an  event  (91057061428.22)  that  6 
occurred  roughly  85  m  to  the  northeast  of  the  array.  Note  the  phase  arriving  2.25 

s  after  the  onset  energy.  The  recordings  have  been  band-passed  between  4  and  6 
Hz. 

5  3-component  image  calculated  using  event  91057061428.22  on  the  vertical  7 
trending  14°  East  of  North.  Energy  is  concentrated  at  a  depth  of  30  km  roughly 

20  to  30  km  to  the  north  of  the  PFBA. 

6  Three  dimensional  image  calculated  using  event  91057061428.22.  Note  the  8 
distribution  of  energy  in  space  along  the  2.25  s  isotime  curve.  The  beaded  nature 

of  the  energy  distribution  is  an  artifact  of  the  gridding  end  with  and  is  not  due  to 
sidelobes  in  the  Impulse  Response. 

7  Landers  and  Superstition  Hills  aftershocks  are  shown  on  the  left  (with  ANZA  13 
network  stations  represented  by  the  triangles).  On  the  right  are  the  two  aftershock 
sequences  in  range  and  depth. 

8  On  the  upper  left  we  show  a  vertical  suite  of  WI  earthquake  synthetics.  14 
Superimposed  on  the  synthetics  are  travel  time  curves  for  pP,  sP,  pPmP  and 
sPmP.  On  the  upper  right  we  show  the  same  traces  after  conversion  to 
envelopes.  On  the  bottom  we  show  the  travel  time  curves  resulting  from  binning 
sources  in  range  from  0.8  to  0.9°.  The  binning  has  the  most  significant  affect  on 

the  PmP  phases. 

9  An  event  column  based  on  Landers  aftershocks  occurring  0.9  to  1.1°  from  the  15 
station.  Superimposed  on  this  figure  are  travel  time  curves  for  pP,  sP,  pPmP  and 
sPmP.  To  the  right  is  a  histrogram  indicating  the  number  of  recordings  in  each 
bin. 

10  This  figure  is  the  same  as  Figure  9  except  it  uses  Landers  aftershocks  located  0.8  16 

to  0.9°  from  the  station. 

1 1  This  figure  is  the  same  as  Figure  9  except  it  uses  Superstition  Hills  aftershocks  17 
located  0.8  to  0.9°  from  the  station. 
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12  All  single-  to  master-trace  cross  correlations  are  shown  here  (for  the  event  18 
columns  shown  in  Figures  9  and  10)  and  are  represented  by  the  triangles.  The 
squares  represent  the  best  single-station  cross-correlations  (i.e.  the  single-station 
depth  solutions).  These  solutions  cluster  around  a  zero  residual. 

13  ANZA  network  depth  residuals  obtained  using  the  event  columns  displayed  in  19 
Figures  9  and  10. 

14  The  NRDC  and  KNET  deployments.  The  NRDC  figure  shows  event  locations.  22 
The  highlighted  events  include  CHI  and  CH2  (the  calibration  explosions)  and 
EVC  (a  quarry  blast). 

15  Two  sonograms  (with  time-series)  calculated  using  NRDC  recordings.  On  the  23 
top  is  the  sonogram  of  CH2  recorded  at  Bayanaul.  On  the  bottom  is  the 
sonogram  of  EVC  also  recorded  at  Bayanaul. 

16  Four  binary  sonograms  calculated  using  NRDC  recordings.  On  the  top  row  are  24 
binary  sonograms  calculated  using  the  recordings  displayed  in  Figure  15.  CH2  is 

on  the  left,  EVC  is  on  the  right.  On  the  bottom  row  are  binary  sonograms  of  EVC 
calculated  using  the  N-S  and  E-W  Bayanaul  recordings. 

17  The  results  of  applying  the  similarity  discriminant  to  the  NRDC  dataset.  In  this  25 
figure  events  1  to  15  are  believed  to  be  quarry  blasts,  16  to  18  are  the  calibration 
explosions.  On  the  left  we  show  single  station  solutions,  the  full  network 
solutions  (using  the  borehole  high-gain  channel)  are  shown  on  the  upper  right 
(with  the  number  of  stations  displayed  below). 

18  The  KNET  and  the  events  used  in  Figures  19  to  21.  The  three  circles  represent  26 
the  three  event  populations.  The  upper  right  circle  encloses  the  limestone  quarry 
location  that  gave  rise  to  the  50  quarry  blasts  used  in  the  calibration  and  testing  of 

the  similarity  discriminant.  In  the  circle  on  the  lower  left  are  the  earthquakes 
(aftershocks  to  the  Mw  7.4  Suusamyr  thrust  earthquake)  used  in  calibrating  the 
discriminant  in  this  area.  The  middle  circle  encloses  the  earthquakes  that  were 
used  to  test  the  discriminant. 

19  Single  station  solutions  and  ranges.  Events  1  to  25  are  believed  to  be  quarry  27 
blasts.  The  remaining  events  are  aftershocks  to  the  Suusamyr  thrust  event. 

20  Network  solutions  using  the  same  (calibration)  events  considered  in  the  preceding  28 
figure.  On  the  bottom  we  display  the  number  of  stations  used  in  each  solution. 

The  earthquake  solutions  follow  a  normal  distribution.  The  horizontal  lines 
represent  the  mean  and  -i-/-  1,2  and  3  standard  deviations. 

21  Network  solutions  using  the  test  events  (see  Figure  18).  The  horizontal  lines  in  29 
the  top  figure  are  reproduced  from  Figure  20  (the  calibration  events). 


Appendix 

A1  Geographic  locations  of  receivers  and  events  in  the  1989  NRDC  dataset.  49 

A2  Seismogram  resulting  from  the  single-chemical  explosion  Chemex  2  recorded  at  a  50 
range  of  157  km  by  the  vertical  component  surface  seismometer  at  Bayanaul  and 
corresponding  sonogram. 

A3  Seismogram  resulting  from  event  c  recorded  at  a  range  of  264  km  by  the  vertical  5 1 
component  surface  seismometer  at  Bayanaul  and  corresponding  sonogram.  This 
event  is  believed  to  be  ripple  fired.  The  sonogram  matrix  contains  a  clear  time- 
independent  spectral  signature. 
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A4  Original  spectral  estimate  (solid)  and  two  filtered  versions,  one  relatively  52 
unsmoothed  (fine  dash)  and  the  other  heavily  smoothed  (coarse  dash).  This 
figure  is  intended  to  illustrate  the  means  by  which  we  reduce  each  spectral 
estimate  to  a  binary  spectral  estimate.  In  the  following  figures  regions  of  locally 
high  power  are  represented  by  a  +1  (white).  Regions  deficient  in  power  are 
represented  by  a  -1  (black). 

A5  Vertical  component  binary  sonogram  matrix  corresponding  to  Chemex  2  recorded  53 
at  Bayanaul  (see  Figure  A2).  The  first  35  seconds  of  coda  after  the  compressional 
onset  are  represented  in  this  and  the  next  binary  sonogram. 

A6  Vertical  component  binary  sonogram  matrix  corresponding  to  event  c  recorded  at  54 
Bayanaul  (see  Figure  A3).  Note  the  clear  time-independent  pattern. 
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Site  Characterization  with  Small-Aperture  Arrays 


1.  INTRODUCTION 

An  important  aspect  of  nuclear  monitoring  is  to  understand  the  generation  of  seismic  coda,  particularly 
the  near-station  mechanisms.  In  some  instances,  secondary  arrivals,  due  to  near  station  processes  (inch 
scattering  from  two-  or  three-dimensional  heterogeneities),  may  be  incorrectly  interpreted  as  direct 
arrivals,  including  so-called  “depth  phases”,  particularly  when  they  are  detected  by  single  stations. 
Studies  quantifying  the  effects  of  scattering  have  mostly  been  statistical  and  have  dealt  with  the  influ¬ 
ence  of  small-scale  random  scatterers  (e.g.  Aki,  1969,  Aki  &  Chouet,  1975).  Some  deterministic  stud¬ 
ies  (Key,  1967;  Key,  1968;  Gupta  et  al.,  1990;  Gupta,  Lynnes  &  Wagner,  1990;  Bannister,  Husebye  & 
Ruud,  1990;  Lay,  1987;  Lynnes  &  Lay,  1989;  Hedlin  et  al,  1994  Revenaugh,  1995)  have  produced 
compelling  evidence  that  large  crustal  inhomogeneities  or  topographic  features  are  capable  of  produc¬ 
ing  significant  amounts  of  scattered  energy,  often  in  the  form  of  large,  identifiable  seismic  phases. 

1.1  Deterministic  Characterizations  via  Kirchoff  Imaging 

Assuming  that  the  velocity  structure  of  the  Earth  volume  of  interest  is  reasonably  well  constrained  and 
that  the  seismic  signals  scattered  from  prominent  heterogeneities  are  recorded  by  a  spatial  deployment 
of  seismometers  it  should  be  possible  to  “reverse”  the  propagation  and  image  the  points  at  which  the 
energy  was  scattered  and  re-directed  toward  the  recording  station(s).  There  are  a  number  of  classic 
imaging  techniques  available  to  “de  propagate”  scattered  energy.  Surface  integration  or  Kirchoff  tech¬ 
niques  are  widely  used  by  Oil  Exploration  firms  to  process  seismic  reflection  data  and  probe  the  sub¬ 
surface  for  hydrocarbon  bearing  structures  or  strata.  In  essence  the  velocity  model  is  used  to  predict 
when,  at  a  given  station,  a  scattered  arrival  should  be  expected  from  a  particular  location  in  the  Earth. 
By  predicting  arrival  times  at  all  recording  stations  and  integrating  the  appropriate  portions  of  the 
records,  any  phases  arriving  from  a  particular  location  in  the  Earth  can  be  constructively  interfered  and 
enhanced.  By  systematically  repeating  this  time-shifting  and  integration  process  for  small  sub-vol¬ 
umes  within  a  volume  of  interest  prominent  scatterers,  or  heterogeneities,  can  be  imaged. 

1.2  Locale  for  the  analysis 

For  this  work  we  have  chosen  seismic  recordings  made  in  a  structurally  complex  area  -  southern  Cali¬ 
fornia.  Although  most  previous  deterministic  analyses  have  focused  on  scattering  at  the  free-surface, 
in  this  work  we  scan  the  entire  crustal  volume  near  the  receiving  station  in  search  of  prominent  scatter¬ 
ers.  Although  multiple  scattering  and  the  chance  coincidence  of  energy  from  small  scale  scatterers  can 
give  rise  to  prominent  coda  phases,  this  study  that  we  will  be  describing  specifically  seeks  the  product 
of  single-discrete  scattering  interactions  at  intermediate  (1  to  10  Hz)  frequencies. 

2  INITIAL  LOOK  AT  DATA 

2.1  ANZA  network  and  PFBA  recordings 

Before  discussing  our  attempts  to  deterministically  associate  specific  coda  phases  with  specific  loca¬ 
tions  in  the  Earth  weTl  conduct  a  brief  review  of  time  series  recorded  in  the  area  of  current  interest.  In 
this  paper  we  are  looking  at  local  recordings  of  seismic  events  in  southern  California  made  by  two 
deployments.  The  ANZA  broadband  network  is  a  10  station  3-component  network  that  utilizes  STS-2 
broadband  sensors  and  24-bit  dataloggers  generating  triggered  (250  sps)  and  continuous  (20  sps)  data 
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Figure  1:  On  the  left  we  show  the  PFBA  (black  circle)  amid  southern  California  to¬ 
pography  and  faults.  The  map  on  the  lower  right  shows  the  PFB  A  and  ANZA  networic 
stations  (black  triangles).  The  PFBA  (upper  right)  is  a  28  element  6  km  aperture  array 
which  was  deployed  for  3  months  in  early  1991.  Each  element  consists  of  3  orthogo¬ 
nal  STS-2  sensors  (.0083  to  40  Hz).  The  Salton  Trough  lies  to  the  north  and  east  of 
these  deployments. 


streams.  The  network  is  located  in  the  Peninsular  Ranges  batholith  in  southern  California  and  was 
emplaced  to  monitor  local  and  regional  seismicity  in  southern  California  (Figure  1).  This  location  has 
proven  to  be  particularly  useful  since  the  batholith  provides  exceptionally  low-loss  and  homogeneous 
transmission  paths  (by  California  standards).  The  network  station  at  the  Pinon  Flat  Observatory  (PFO) 
is  particularly  useful  since  it  is  deployed  at  the  site  that  was  used  for  a  temporary  array,  the  Pinon  Flat 
Broadband  Array  (PFBA).  The  PFBA  was  deployed  during  the  winter  months  of  1991  and  consisted 
of  28  stations  located  in  concentric  rings  with  a  maximum  aperture  of  6  km  (Figure  1).  Each  station  in 
the  array  was  comprised  of  3  orthogonal  broadband  (0.0083  to  40  Hz)  sensors.  During  its  short  life  the 
array  made  triggered  (100  sa/s)  recordings  of  130  local  and  140  regional  and  teleseismic  events  as  well 
as  continuous  (20  sa/s)  recordings  of  these,  and  other  smaller,  events. 
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Figure  2:  Record  section  of  events  recorded  by  the  vertical  component  of  the  ANZA  network  station 
PFO  (located  at  the  intersection  of  the  crosshairs  in  the  map).  The  events  occurred  between  0.71  and 
0.78°  from  the  station.  This  record  section  includes  events  at  depths  ranging  from  0.8  to  5  km,  Ml 
ranges  from  1.9  to  3.1.  Despite  the  similar  propagation  lengths  these  seismic  records  display  little  in 
common. 


For  an  initial  look  at  waveforms  we  display  (in  Figure  2)  vertical  component  recordings  of  local  events 
made  by  the  network  station  at  PFO.  This  figure  contains  events  that  occurred  between  0.71  and  0.78° 
from  PFO  at  depths  ranging  from  .8  to  5  km  and  had  local  magnitudes  of  1.9  to  3.1.  The  time  series, 
which  are  aligned  on  the  first  arrivals,  are  remarkably  dissimilar.  Some  records  (for  example  those  at 
a  back  azimuth  from  PFO  of  35°)  contain  a  discrete  high-amplitude  arrival  roughly  2  s  after  onset. 
Other  events  contain  a  similar  arrival  elsewhere  in  the  P  wave  coda  {e.g.  earlier  at  a  back  azimuth  of 
75°).  Some  events  contain  no  distinct  coda  phases.  Some  differences  could  be  due  to  variations  in  the 
Moho  depth,  differences  in  event  range  from  PFO,  depth  or  focal  mechanism.  Examination  of  this 
figure  alone  might  lead  one  to  conclude  that  the  high  amplitude  coda  phases  are  simply  near-critical 
Moho  reflections;  however,  a  closer  look  at  the  data  suggests  otherwise. 

The  ANZA  network  made  triggered  and  continuous  mode  recordings  of  the  Landers  event  and  158 
aftershocks  above  M  4.0  (Figure  3).  This  linear  aftershock  trend  gives  us  an  opportunity  to  view  the 
behavior  of  coda  phases  with  range  since  the  ANZA  network  stations  recorded  the  events  at  ranges 
from  0.2  to  1.2°.  It  is  likely  that  the  aftershocks  have  similar  focal  mechanisms.  In  Figure  3  we  show 
record  sections  of  aftershocks  ranging  from  0.6  to  0.9°  from  three  stations  in  the  network  (KNW,  FRD 
and  PFO).  Although  there  is  some  indication  of  a  late  phase  in  the  KNW  and  FRD  record  sections  only 
the  PFO  records  consistently  contain  a  bright  coda  phase.  This  phase  appears  to  “turn  on”  at  roughly 
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Figure  3  In  this  figure  we  display  vertical  component  recordings  of  Landers  aftershocks 
(above  Ml  3.5)  recorded  by  three  of  the  ANZA  network  stations.  The  Pinon  Flat  station 
shows  that  an  obvious  phase  (delayed  2  s  after  the  onset)  appears  shortly  before  .7°  and 
suggests  that  it  disappears  beyond  0.85°.  The  sudden  appearance  of  this  phase  suggests  that 
it  is  a  bounce  from  the  Moho  that  reaches  the  critical  angle  by  0.7°.  The  other  two  stations 
suggest  that  the  energy  is  not  due  to  a  simple  bounce  since  there  is  little  evidence  of  a  late 
phase  at  these  locations. 
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0.7°  and  might  not  be  present  beyond  0.85°.  Not  only  is  this  phase  capricious  among  events  distributed 
azimuthally  but  also  among  events  distributed  in  range.  Furthermore,  depending  on  where  an  event  is 
recorded  the  phase  might,  or  might  not,  be  present. 

Although  the  PFBA  was  deployed  after  the  Landers  aftershock  sequence  this  array  recorded  some  of 
the  events  considered  in  Figure  2,  including  the  two  mentioned  above  located  at  a  back  azimuth  of  35°. 
An/-k  analysis  of  PFBA  recordings  of  events  located  to  the  northeast  of  the  array/network  (summa¬ 
rized  in  Hedlin  &  Vemon,  1995)  show  that  the  onset  phase  crosses  the  array  from  the  azimuth  of  the 
source  at  a  phase  velocity  of  6  to  7  km/s  while  the  later  phase  crosses  with  an  approach  angle  biased  to 
the  north  at  a  phase  velocity  of  10  to  11  km/s. 

This  preliminary  look  at  data  has  given  us  some  idea  of  the  extreme  complexity  of  local  coda  recorded 
in  southern  California  and  has  led  to  several  questions:  1)  Why  do  some  local  coda  contain  bright 
phases  while  others,  from  nearby  events,  do  not?  2)  Why  do  certain  stations  seem  to  receive  bright  coda 
phases  while  others,  nearby,  do  not?  3)  Why  do  these  late  phases  sometimes  arrived  from  biased  direc¬ 
tions?  To  further  our  understanding  of  these  coda  phases  we  will  use  the  3-component  PFBA  to  scan 
the  crust  and  upper  mantle  and  image  the  point(s)  at  which  they  are  redirected  toward  the  receiver. 

3.  BRIEF  REVIEW  OF  THE  IMAGING  TECHNIQUE 

In  calculating  images  of  crustal  scatterers  we  follow  a  simple  hypothesis.  We  invoke  the  first  Born 
approximation  and  thus  assume  that  coda  energy  is  solely  due  to  single  scattering  at  locations  within 
the  Earth  or  at  the  surface.  By  following  this  assumption  in  our  processing  we  suppress  contributions 
made  to  coda  by  higher  order  scattering.  We  are  deliberately  attempting  to  deterministically  character¬ 
ize  single  scattering  at  the  more  prominent  crustal  scatterers  or  interfaces.  Initial  versions  of  this 
technique  have  been  described  in  earlier  papers  (Hedlin,  1991;  Hedlin  et  al.  1991,  1994)  and  so  only  a 
brief  review  and  a  summary  of  the  improvements  will  be  given  here.  In  our  analysis  of  free-surface 
scatterers  in  southern  Norway  we  used  single  (vertical)  component  station  recordings  of  teleseisms 
made  by  the  small-aperture  NORESS  array.  Imaging  the  region  surrounding  the  array  was  accom¬ 
plished  simply  by  conducting  a  systemmatic  (grid)  search  for  scatterers.  At  each  point  we  would 
beamform  for  scattered  waves  coming  from  the  candidate  scatterer  and  use  simple  travel  time  calcula¬ 
tions  to  predict  where  in  the  array  records  any  scattered  waves  should  appear.  In  essence  we  assume 
that  a  scatterer  exists  at  a  location  in  a  grid  and  then  we  stack  the  records  to  constructively  interfere 
single-scattered  signals.  This  approach  is  a  time-surface  integration,  or  Kirchoff  migration,  technique 
(Yilmaz,  1987).  Hedlin  et  al.  (1991;  1994)  used  a  similar  approach  to  image  P  to  Rg  scatterers  in 
Norway  and  found  evidence  that  energy  was  being  scattered  at  sharp  topographic  relief. 

In  this  paper  we  extend  our  earlier  work  in  several  ways.  Instead  of  restricting  our  analysis  to  surface 
scatterers  excited  by  teleseismic  wavetrains  we  have  adapted  the  technique  to  accommodate  scatterers 
located  anywhere  in  a  layered  earth  excited  by  local,  regional  or  teleseismic  events.  We  have  intro¬ 
duced  a  polarization  analysis  into  the  technique  that  takes  advantage  of  3-C  array  data.  To  backproject 
the  energy  into  a  volume  of  the  crust  we  look  at  individual  sub-volumes  in  sequence.  Raytracing 
through  an  assumed  Earth  model  associates  each  sub-volume  with  a  propagation  vector  and  with  the 
appropriate  time  delay  that  should  exist  between  the  parent  and  a  single-scattered  phase.  The  3-com¬ 
ponent  array  data  are  then  tuned  to  this  vector  by  steering  each  3-component  station  for  the  direction  of 
particle  motion  of  interest  (e.g.  longitudinal,  transverse)  and  then  through  beamforming  these  steered 
channels  are  summed  together  to  yield  a  single  trace  dedicated  to  energy  propagating  from  the  subvolume. 
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Figure  4:  Vertical  component  triggered-mode  recordings  of  an  event  (91057061428.22)  that 
occurred  roughly  85  m  to  the  northeast  of  the  array.  Note  the  phase  arriving  2.25  s  after  the  onset 
energy.  The  recordings  have  been  band-passed  between  4  and  6  Hz. 

The  time-delay  calculated  for  the  sub-volume  determines  the  portion  of  the  tuned  time-series  that  is 
attributed  to  it.  This  process,  which  is  akin  to  Kirchoff  migration,  yields  an  image  of  coherent,  first- 
order  scatterers  that  exist  within  the  Earth  or  at  the  free  surface.  Having  3-component  data  obviates  the 
need  for  an  obliquity  factor  that  is  used  when  Kirchoff  migration  is  applied  to  vertical  component 
reflection  data  (Yilmaz,  1987).  To  date  we  have  restricted  the  analysis  to  1 -Dimensional  Earth  models 
and  apply  the  appropriate  geometrical  spreading  factor. 

4.  DATA  ANALYSIS 

4.1  Imaging  analysis  of  PFBA  data 

To  demonstrate  this  technique  weTl  consider  the  bright  coda  phase  discussed  in  section  2.  We  have 
analyzed  triggered  and  continuous  mode  recordings  of  local  and  regional  events  that  occurred  to  the 
northeast  of  the  array.  One  event  (which  occurred  on  day  057)  occurred  85  km  from  the  array  at  a  back- 
azimuth  of  35°.  The  P  and  S  onsets  occur  at  roughly  10  and  20  s  respectively.  The  coda  contains  a 
significant  phase  delayed  2.25  s  after  the  P  wave  onset  (Figure  4).  As  discussed  in  section  2  f-k 
analysis  indicates  that  the  onset  phase  arrives  from  the  azimuth  of  the  event,  at  a  phase  velocity  of  6.2 
km/s  while  the  later  phase  arrives  from  a  back  azimuth  of  14°  at  12  km/s.  In  Figure  5  we  show  a  profile 
through  the  crust  that  passes  through  the  array  and  bears  14°  East  of  North.  In  the  upper  figure  we 
display  the  time  delays  we  predict  should  exist  between  the  arrivals  of  the  energy  directly  from  the 
event  and  via  single  P-P  scattering  interactions  that  occur  in  this  vertical  plane.  Given  that  our  assump¬ 
tions  (regarding  the  velocity  model  and  the  scattering)  are  reasonable  the  source  of  the  late  energy 
should  lie  somewhere  on  the  2.25  s  curve  (dashed).  The  middle  third  of  Figure  5  shows  the  slowness 
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Figure  6:  Three  dimensional  image  calculated  using  event  91057061428.22.  Note  the  distribu¬ 
tion  of  energy  in  space  along  the  2.25  s  iso  time  curve.  The  beaded  nature  of  the  energy  distiibu- 
tion  is  an  artifact  of  the  gridding  end  with  and  is  not  due  to  sidelobes  in  the  Impulse  Response. 

we  expect  energy  resulting  from  P-P  scattering  in  this  plane  should  have.  This  figure  allows  us  to 
predict  where  in  the  plane  the  scatterer  should  lie  -  at  the  intersection  of  the  dashed  curves  in  the  upper 
and  lower  figures.  At  the  bottom  of  Figure  5  we  display  an  image  of  scattering  intensity  calculated 
from  this  event  using  the  19  available  3-C  stations.  The  energy  lies  close  to  the  2.25  s  curve  and  despite 
the  imperfect  slowness-resolution  it  is  clear  that  an  isolated  source  exists  at  roughly  30  km  depth  and 
roughly  20  km  to  the  NNE  of  the  array.  In  Figure  6  we  display  a  3-D  image  of  the  crust  beneath  the 
array  obtained  with  our  3C  migration  technique.  The  image  extends  from  21  km  south  and  west  of  the 
array  to  62  km  north  and  east  between  a  depth  of  20  to  40  km  and  is  constructed  from  41  x  41  x  31  grid 
points  (in  x,y  and  z).  The  comer  has  been  removed  to  reveal  the  concentration  of  energy  at  30  km 
depth.  In  this  figure  the  effect  of  the  limited  resolution  is  indicated  by  the  smearing  of  energy  along  the 
ellipsoidal  surface.  The  beading  is  an  artifact  of  gridding  associated  with  the  calculation  of  the  image. 
The  proximity  of  this  feature  to  the  San  Andreas  fault  is  intriguing  (it  is  just  south  of  the  point  of 
intersection  with  the  Landers  fault  zone).  The  depth  is  puzzling,  however,  since  at  30  km  it  is  well 
below  the  seismogenic  zone. 

5.  DISCUSSION 

Since  the  rays  involved  cross  a  complicated  crust  (which  includes  the  Coachella  valley)  there  are  a 
number  of  possible  explanations.  For  example  the  late  arrival  might  be  a  deeply  penetrating,  normally 
refracted,  ray  that  has  passed  under  the  trough.  It  is  possible  that  the  late  phase  is  not  due  to  normal 
refraction  but  has  been  caused  by  scattering  from  a  deeply  buried  feature  at  the  northwestern  end  of  the 
valley.  Considering  that  a  scattering  process  is  unlikely  to  yield  such  a  high  amplitude  arrival  it  seems 
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most  likely  that  this  phase  has  resulted  from  a  bounce  from  a  continuous  interface.  At  a  depth  of  30  km 
this  feature  lies  at  the  depth  of  the  Moho  {Hadley  &  Kanamori,  1977).  The  amplitude  is  suggestive  of 
a  near-critical  angle  reflection  from  the  Moho,  it  is  hard  to  imagine  any  scattering  process  that  will 
yield  that  much  energy.  The  fact  that  the  phase  is  often  absent,  and  when  present  often  arrives  from  a 
biased  direction,  requires  that  the  phase  is  not  due  to  a  simple  bounce  from  a  flat  Moho.  Although  the 
velocity  model  we  have  used  for  the  raytracing  is  1 -dimensional  there  is  clearly  a  third  dimension  to 
the  velocity  structure  and  it  seems  most  likely  that  the  phase  is  bouncing  from  a  highly  contorted 
Moho.  Although  PFO  lies  on  a  stable,  Mesozoic  age,  granitic  batholith  the  larger  region  is  tectonically 
active  -  as  evidenced  by  the  topography  displayed  in  Figure  1 .  The  Moho  discontinuity  does  not  likely 
approximate  a  horizontal  plane  but  is  likely  contorted.  The  feature  we  have  identified  lies  just  south  of 
the  San  Andreas  fault  zone  at  the  northern  end  of  the  Coachella  Valley.  Although  no  study  has  clearly 
indicated  the  orientation  of  the  Moho  surface  in  this  area  it  is  widely  believed  (e.g.  Sung  &  Jackson, 
1992;  Zhao,  Kanamori  &  Humphrey,  1994;  Baker  et  al,  1995)  that  the  Moho  tilts  down  to  the  north¬ 
west  in  this  area  as  leaves  the  Salton  trough.  The  regional  tilt  inferred  from  these  studies  is  inconsistent 
with  the  tilt  necessary  for  Moho  bounces  to  be  directed  toward  PFO  from  events  to  the  northeast. 
Jachens  &  Griscom  ( 1985)  mapped  isostatic  gravity  residuals  assuming  Airy  compensation  at  the  Moho 
and  found  no  significant  residuals.  This  indicates  that  crustal  thickness  can  be  inferred  from  free- 
surface  topography  -  i.e.  an  inverse  correlation  should  exist  between  topography  on  the  Moho  and  at 
the  free-surface.  As  shown  in  Figure  1  the  topography  drops  sharply,  from  northeast  to  southwest, 
across  the  San  Andreas  fault  north  of  the  ANZA  network  as  the  transverse  ranges  contact  the  Salton 
Trough.  The  topography  rises  sharply  just  south  of  the  trough  in  the  vicinity  of  the  network.  If  this 
topography  is  mirrored  by  topography  on  the  Moho  along  the  same  profile  the  Moho  would  rise  sharply 
at  the  San  Andreas  fault  (assuming  a  density  contrast  at  the  Moho  of  0.6  g/cc  and  a  crustal  density  of 
2.7  g/cc  a  1000  m  drop  at  the  free-surface  would  correspond  to  a  4500  m  rise  in  the  Moho).  A  sharp  rise 
in  the  Moho  trending  northwest-southeast  is  consistent  with  the  topography  that  would  be  necessary  to 
reflect  energy  arriving  from  the  northeast  toward  PFO. 

6.  CONCLUSIONS 

This  project  was  initiated  with  the  goal  of  identifying  prominent,  buried,  scatterers  that  give  rise  to 
identifiable  coda  phases  that  might  be  confused  with  near  source  phases  {e.g.  those  used  to  constrain 
depth).  We  have  observed  extreme  variability  in  seismic  coda  recorded  at  local  and  regional  ranges  by 
ANZA  and  the  PFBAin  southern  California.  We  believe  that  the  variability  is  not  solely  due  to  scatter¬ 
ing  from  volume  heterogeneities  but  to  the  interaction  of  seismic  energy  with  Moho  topography.  This 
work  suggests  that  any  analysis  of  event  location  that  involves  interface  bounce  phases  (particularly 
those  occurring  in  tectonically  active  areas)  should  consider  the  possibility  of  bias  in  the  arrival  direc¬ 
tion  and  arrival  time.  Although  this  study  has  been  used  to  describe  and  demonstrate  the  imaging 
technique  these  results  do  not  represent  a  comprehensive  search  for  heterogeneities  in  the  vicinity  of 
Pinon  Flat. 
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Empirical  Characterization  of  Regional  Seismic  Coda  -  Constraining  Depth 


1.  INTRODUCTION 

We  use  network  recordings  of  local  aftershocks  to  empirically  characterize  local  crustal  wave  propaga¬ 
tion  and  look  for  robust  coda  phases  that  are  dependent  on  source  depth.  These  empirical  results  form 
the  basis  for  a  technique  for  constraining  source  depth  in  seismically  active  areas. 

To  search  for  robust  depth  phases  we  group  aftershocks  into  bins  dependent  on  “known”  event  depth 
and  the  range  to  the  recording  station.  Considering  each  bin  individually  we  stack  the  recordings  after 
converting  each  trace  to  an  envelope  and  aligning  all  traces  on  the  first  arrival.  By  repeating  this 
procedure  for  a  vertical  set  of  bins  we  construct  an  “event  column”.  Using  ANZA  network  recordings 
of  Landers  and  Superstition  Hills  aftershock  sequences  we  have  constructed  vertical  component  event 
columns  at  ranges  from  0.6  to  l.T  spanning  depths  from  the  free  surface  to  10  km.  In  the  resulting 
displays,  which  show  the  dependence  of  de-phased  coda  on  range  and  depth,  we  see  clear  evidence  of 
PmP  (at  ranges  from  0.7  to  1.1°)  but  no  evidence  of  any  free-surface  bounces  at  any  of  these  ranges. 
Although  synthetics  suggest  that,  at  the  ranges  considered,  pPg  and  pPmP  should  have  very  small 
amplitude  (relative  to  Pg  and  PmP)  it  is  surprising  that  no  robust  sPg  and  sPmP  phases  are  evident. 

In  spite  of  the  presence  of  only  a  single  depth  phase  (which  does  not  exhibit  a  rapid  moveout  with 
respect  to  Pg)  we  have  attempted  to  use  these  “event  columns”  to  constrain  event  depth.  A  single¬ 
station  event  depth  is  estimated  by  cross-correlating  the  recording  of  the  event  with  the  appropriate 
event  column  (after  removing  the  trace  being  considered  from  the  event  column).  The  network  depth 
estimate  is  obtained  from  a  weighted  average  of  all  single-station  depth  estimates  (with  preference 
being  given  to  those  stations  with  the  highest  single-trace  to  master-trace  cross-correlations).  From  an 
analysis  of  43  Landers  aftershocks  we  have  found  that  the  cross-correlations  do  not  depend  strongly  on 
depth.  The  depth  residuals  (“known”  depth  -  our  depth  estimate)  have  zero  mean  with  a  standard 
deviation  of  2.5  km.  This  residual  decreases  as  the  number  of  stations  used  increases.  With  7  stations 
“voting”  the  residuals  have  a  standard  deviation  of  2.0  km.  In  areas  where  free-surface  depth  phases 
are  robust  we  expect  that  this  approach  would  yield  smaller  depth  residuals. 

2.  SYNTHETIC  AND  RECORDED  DATA 

2.1  Wavenumber  Integration 

In  this  section  we  use  Wavenumber  Integration  (WI  -  Apsel  and  Luco,  1983;  Luco  and  Apsel,  1983)  to 
synthesize  local  and  regional  waveforms. 

2.2  ANZA  network  recordings 

We  have  examined  local  recordings  of  seismic  events  in  southern  California  made  by  the  ANZA  net¬ 
work.  The  ANZA  network  made  triggered  and  continuous  mode  recordings  of  the  June  28  1992  Mw 

7.3  Landers  event  and  158  aftershocks  above  M  4.0  (Figure  7).  This  linear  aftershock  trend  gives  us  an 
opportunity  to  view  the  behavior  of  coda  phases  with  range  since  the  ANZA  network  stations  recorded 
the  events  at  ranges  from  0.28  to  1 .28°.  The  Mw  6.0  Superstition  Hills  main  shock  occurred  in  Novem¬ 
ber,  1987  and  was  followed  by  35  events  above  Ml  3.5.  The  vast  majority  of  these  events  occurred 
within  10  km  of  the  free-surface. 
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Figure  7:  Landers  and  Superstition  Hills  aftershocks  are  shown  on  the  left  (with  ANZA 
network  stations  represented  by  the  triangles).  On  the  right  are  the  two  aftershock  se¬ 
quences  in  range  and  depth.. 


3.  PRELIMINARY  ANALYSIS 

Prior  to  an  analysis  of  recorded  data  we  have  considered  used  WI  to  predict  the  response  of  the  south¬ 
ern  California  crustal  model  {Zhao  &  Helmberger,  1994)  to  various  dislocation  and  dilatational  sources. 
We  have  calculated  a  suite  of  synthetics  using  WI  for  dislocation  sources  ranging  from  0.5  to  9.5  km 
deep  at  a  range  of  1°  (Figure  8).  The  focal  mechanism  used  was  that  of  a  vertical  right-lateral  strike  slip 
dislocation  on  trend  with  the  Landers  aftershock  sequence  recorded  at  the  site  of  station  PFO  in  ANZA. 
In  this  figure  we  display  travel  times  for  Pg,  pPg,  sPg,  pPmP  and  sPmP  and  find  that  the  pPg  and 
pPmP  phases  are  not  discernible.  At  this  range  we  expect  little  energy  in  the  “little  P”  phases  (due  to 
inefficient  conversion  of  upgoing  compressional  energy  at  the  free  surface).  These  dislocation  sources 
generate  enough  shear  wave  energy  to  yield  significant  sPg  and  sPmP  -  these  phases  are  comparable  in 
amplitude  to  Pg  and  PmP  -  in  addition  to  a  number  of  higher  order  crustal  resonance  phases.  The  “little 
S”  phases  are  potentially  very  useful  for  depth  determination  because  they  move  out  rapidly  with 
respect  to  Pg. 

4.  EVENT  COLUMNS 

By  binning  and  stacking  event  recordings  (after  aligning  on  the  P  onset  and  converting  the  traces  to 
envelopes)  we  can  empirically  characterize  crustal  wave  propagation  and  look  for  coda  variations  with 
event  depth  and  source-receiver  range.  We  expect  some  loss  in  temporal  resolution  will  result  from 
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Figure  8:  On  the  upper  left  we  show  a  vertical  suite  of  WI  earthquake  synthetics.  Super¬ 
imposed  on  the  synthetics  are  travel  time  curves  for  pP,  sP,  pPmP  and  sPmP.  On  the 
upper  right  we  show  the  same  traces  after  conversion  to  envelopes.  On  the  bottom  we 
show  the  travel  time  curves  resulting  from  binning  sources  in  range  from  0.8  to  0.9°. 
The  binning  has  the  most  significant  affect  on  the  PmP  phases. 
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Figure  9;  An  event  column  based  on  Landers  aftershocks  occurring  0.9  to 
1.1°  from  the  station.  Superimposed  on  this  figure  are  travel  time  curves 
for  pP,  sP,  pPmP  and  sPmP.  To  the  right  is  a  histrogram  indicating  the 
number  of  recordings  in  each  bin. 

converting  to  envelopes  and  binning  (Figure  8).  The  binning  should  have  the  greatest  effect  on  the 
PmP  phases. 

Actual  event  columns  are  shown  in  Figures  9  to  1 1 .  The  first  and  second  were  calculated  using  Landers 
aftershocks  from  0.9  to  1.1°  and  from  0.8  to  0.9°.  The  third  column  used  aftershocks  to  the  Supersti¬ 
tion  Hills  mainshock.  Each  column  is  paired  with  a  histogram  showing  the  number  of  events  that  were 
used  in  each  bin.  The  concentration  of  events  at  6  km  in  Figures  19  and  20  is  not  physical  but  an 
artifact  of  the  software  package  used  to  make  the  original  depth  estimates.  These  events  have  been 
excluded  from  further  consideration. 
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Figure  10:  This  figure  is  the  same  as  Figure  9  except  it  uses  Landers  after¬ 
shocks  located  0.8  to  0.9°  from  the  station. 


We  clearly  see  a  PmP  phase  in  all  the  event  columns.  In  the  Superstition  Hills  event  column  this  phase 
arrives  significantly  earlier  and  suggests  a  thinner  crust  (30  vs  35  km).  There  is  no  evidence  in  any  of 
these  event  columns  of  any  free-surface  bounce  phases  (inch  pP  as  expected,  and  sP  phases). 

5.  DEPTH  ESTIMATION 

In  spite  of  the  absence  of  numerous  robust  coda  depth  phases,  and  the  fact  that  the  one  that  is  present, 
PmP,  shows  little  move-out  with  respect  to  Pg,  we  have  attempted  to  use  these  event  columns  to  esti- 
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Figure  11:  This  figure  is  the  same  as  Figure  9  except  it  uses  Supersti¬ 
tion  Hills  aftershocks  located  0.8  to  0.9°  from  the  station. 


mate  the  depths  of  these  events.  The  scheme  we  have  adopted  is  simply  based  on  the  cross  correlation 
between  a  single-station  recording  of  an  event  and  all  master-traces  in  the  event  column  at  the  source- 
receiver  range.  In  principle,  given  a  vertically  varying  suite  of  master-traces,  this  cross  correlation 
should  be  greatest  at  the  actual  event  depth.  One  advantage  of  this  approach  is  that  it  depends  on 
phases  being  present  in  coda  but  it  is  not  necessary  that  these  phases  are  identified  or  modeled  -  it  is 
simply  based  on  the  shape  of  the  envelope  of  the  entire  waveform.  It  is  assumed  that  the  event  range  is 
constrained  well  enough  to  allow  comparison  with  the  correct  event  column. 

As  an  example  of  this  approach  we  show,  in  Figure  12,  single-  to  master-trace  cross  correlations  as  a 
function  of  the  residual  between  the  known  event  depth  and  depth  of  the  master-trace  (known  -  our 
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Figure  12:  All  single-  to  master-trace  cross  correlations  are  shown  here  (for  the  event  columns 
shown  in  Figures  9  and  10)  and  are  represented  by  the  triangles.  The  squares  represent  the 
best  single-station  cross-correlations  (i.e.  the  single-station  depth  solutions).  These  solutions 
cluster  around  a  zero  residual. 


depth)  for  events  located  in  the  Landers  aftershock  sequence  between  0.8  and  0.9°.  The  single-station 
solutions  (the  depths  at  which  the  cross-correlation  was  greatest)  are  represented  by  the  boxes  and  are 
clearly  clustered  about  a  zero  residual.  The  mean  of  the  solutions  is  -0.89  km  with  a  standard  deviation 
of  3.3  km.  The  standard  deviation  of  the  residuals  is  directly  dependent  on  the  rate  at  which  the  cross¬ 
correlations  change  with  depth.  Due  to  the  absence  of  any  coda  phases  with  rapid  moveout  with- 
respect-to  Pg  this  variation  is  subtle.  Considering  Landers  aftershocks  from  0.9  to  1.1°  from  the  station 
(Figure  12  right  side)  the  result  is  essentially  the  same. 

By  considering  all  network  recordings  of  an  event  simultaneously  the  depth  residuals  can  be  decreased. 
By  simply  forming  a  weighted  average  depth  estimate  (with  weights  given  by  the  single-  to  master- 
trace  cross  correlations)  we  find  that  the  depth  residuals  decrease  with  increasing  numbers  of  stations 
(Figure  13).  With  7  stations  voting  the  standard  deviation  is  2.0  km.  With  two  stations  reporting  the 
residual  is  3.46  km.  The  standard  deviation  for  all  events  located  with  multiple  stations  is  2.5  km. 


6.  CONCLUSIONS 

An  empirical  technique  for  characterizing  local  and  regional  event  coda,  like  the  one  described  above, 
might  be  useful  in  seismically  active  areas  where  the  possibility  of  clandestine  testing  exits.  It  might  be 
sufficient  to  characterize  the  coda  using  a  temporary  deployment,  and  then  using  a  smaller  scale,  more 
permanent,  local  or  regional  network  to  monitor  seismicity  and  contribute  time-series  for  comparison 
to  the  master-traces.  Unless  monitoring  occurs  in  an  area  in  which  robust  free-surface  phases  are 
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Figure  13:  ANZA  network  depth  residuals  obtained  us¬ 
ing  the  event  columns  displayed  in  Figures  9  and  10. 

present  then  a  significant  (~10)  number  of  stations  might  be  necessary  to  yield  the  necessary  depth 
resolution  to  distinguish  deep  earthquakes  from  events  shallow  enough  to  possibly  be  nuclear  tests. 
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A  new  time-frequency  discriminant  based  on  similarity 


1.  INTRODUCTION 

Under  previous  Air  Force  contracts  we  developed  a  discriminant  that  seeks  long-lived  spectral  modu¬ 
lations  in  the  major  phases  and  coda.  Long  lived  modulations  can  be  produced  by  seismic  resonance 
and  by  ripple  firing  and  hence,  given  a  non-resonant  crust,  can  be  diagnostic  of  ripple-fired  events.  In 
order  to  aid  automated  recognition  Hedlin  et  al  { 1989)  developed  a  procedure  in  which  a  binary  sonogram 
is  derived  from  the  original,  spectral,  sonogram  by  the  application  of  filters  which  replace  local  spec¬ 
tral  information  with  a  binary  code  which  simply  reflects  local  spectral  highs  and  lows.  The  binary 
conversion  is  achieved  by  convolving  each  spectral  estimate  within  the  sonogram  with  two  boxcar 
functions  to  yield  two  smoothed  versions  of  the  spectrum  -  one  heavily  smoothed  and  the  other  rela¬ 
tively  unsmoothed.  The  binary  sonogram  simply  results  from  the  difference  between  the  two  sonograms 
with  local  highs  and  lows  represented  by  +/-  1  respectively. 

In  Figure  15  we  show  sonograms  calculated  from  regional  observations  of  a  single-explosion  (top)  and 
a  quarry  blast  (bottom)  recorded  by  the  NRDC  network  in  Kazakhstan  (Figure  14).  The  corresponding 
binary  sonograms  are  shown  in  Figure  16  (top  row).  On  the  bottom  row  we  show  binary  sonograms 
calculated  from  the  North-South  and  East- West  recordings  of  the  quarry  blast.  These  binary  sonograms 
show  a  remarkable  independence  of  the  energy  from  the  recording  component  and  from  time.  Previ¬ 
ously  we  examined  the  time-independence  as  a  discriminant  {Hedlin  et  al,  1989;  Hedlin  et  al,  1990). 
In  this  section  we  will  determine  if  the  independence  from  the  recording  component  might  also  be 
useful  as  a  discriminant. 

2.  DATASETS 

For  this  section  we  have  chosen  to  analyze  two  datasets  collected  in  central  Asia.  The  NRDC  network 
(Figure  14)  was  deployed  in  1987  in  central  Kazakhstan  and  collected  triggered  250  sps  recordings  at 
the  surface  and  in  boreholes  of  three  calibration  explosions,  numerous  local  and  regional  quarry  blasts 
and  a  large  number  of  far-regional/teleseismic  earthquakes. 

The  10-station  telemetered  KNET  was  deployed  in  Kyrgyzstan  in  1991  (Figure  14).  Located  within  an 
aperture  of  200  km  between  the  Kazakh  platform  to  the  north  and  the  Tien  Shan  to  the  south,  this 
network  is  ideal  for  recording  and  studying  events  along  the  seismically  active  southern  flank  of  Eurasia. 
The  network  has  made  100  sps  triggered  and  20  sps  continuous  mode  recordings  of  quarry  blasts  and 
earthquakes  and  local  and  regional  distances  in  addition  to  recording  a  large  number  of  teleseisms. 

2.1  A  priori  event  ID: 

Given  that  we  have  not  had  access  to  “ground-truth”  information,  which  would  positively  identify  the 
type  of  each  event  analyzed,  we  have  relied  on  various  non-seismic  sources  of  information.  Identifica¬ 
tion  of  the  non-calibration  events  in  the  seismically  quiet  Kazakh  platform  were  provided  by  Cliff 
Thurber  who  used  SPOT  photos  to  confirm  mining  activity  at  event  locations.  We  have  found  no 
descriptions  of  shooting  technique  at  these  locations.  Event  identification  in  Kyrgyzstan  was  some¬ 
what  more  difficult.  With  a  hand-held  GPS  unit  Rob  Mellors  identified  an  area  northeast  of  station 
CHM  which  contained  a  large  limestone  quarry  (43.028  N,  74.888  E)  with  shots  typically  occurring  in 
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Figure  14:  The  NRDC  and  KNET  deployments.  The  NRDC  figure  shows  event 
locations.  The  highlighted  events  include  CHI  and  CH2  (the  calibration  explosions) 
and  EVC  (a  quarry  blast). 


Figure  15:  Two  sonograms  (with  time-series)  calculated  using  NRDC  recordings. 
On  the  top  is  the  sonogram  of  CH2  recorded  at  Bayanaul.  On  the  bottom  is  the 
sonogram  of  EVC  also  recorded  at  Bayanaul. 
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Figure  16:  Four  binary  sonograms  calculated  using  NRDC  recordings.  On  the  top  row  are 
binary  sonograms  calculated  using  the  recordings  displayed  in  Figure  15.  CH2  is  on  the  left, 
EVC  is  on  the  right.  On  the  bottom  row  are  binary  sonograms  of  EVC  calculated  using  the  N-S 
and  E-W  Bayanaul  recordings. 


the  late  afternoon.  Recordings  of  some  of  these  events  made  at  CHM  contain  obvious  air  blasts  roughly 
37  s  after  the  P  onset. 

3.  ANALYSIS  OF  KNET  AND  NRDC  DATA 

As  we  noted  in  section  1  ripple  fired  events  are  capable  of  generating  coda  that  is  remarkably  indepen¬ 
dent  of  the  recording  direction.  In  this  section  we  attempt  to  determine  if  this  trait  can  be  used  to 
distinguish  these  events  from  earthquakes  and  explosions.  Using  KNET  and  NRDC  we  consider  dis¬ 
crimination  with  single  stations  and  full  networks. 

The  procedure  that  we  have  developed  to  quantify  the  independence  of  the  time-frequency  patterns 
from  the  recording  direction  is  based  on  a  zero-lag  cross-correlation.  In  Figures  17  to  21  we  display  the 
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Figure  17:  The  results  of  applying  the  similarity  discriminant  to  the  NRDC  dataset.  In  this  figure 
events  1  to  15  are  believed  to  be  quarry  blasts,  16  to  18  are  the  calibration  explosions.  On  the  left 
we  show  single  station  solutions,  the  full  network  solutions  (using  the  borehole  high-gain  chan¬ 
nel)  are  shown  on  the  upper  right  (with  the  number  of  stations  displayed  below). 


average  of  a  three-way  cross-correlation  between  each  pair  of  recording  channels  (i.e.  Z  v  E-W/  Z  v  N- 
S  &  N-S  V  E-W). 

3.1  NRDC: 

Considering  single-station  solutions  the  borehole  sensors  yield  a  better  separation,  the  surface  station 
at  KSU  is  resonant  (Figure  17)  and  thus  has  a  high  degree  of  similarity  between  the  recording  direc¬ 
tions.  Considering  borehole  recordings  the  calibration  explosions  are  tightly  grouped.  There  are  clearly 
two  outliers  among  the  quarry  blasts  (events  1  and  10)  which  are  discussed  below.  The  network  solu- 
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Figure  18:  The  KNET  and  the  events  used  in  Figures  19  to  21.  The  three  circles  represent  the 
three  event  populations.  The  upper  right  circle  encloses  the  limestone  quarry  location  that  gave 
rise  to  the  50  quarry  blasts  used  in  the  calibration  and  testing  of  the  similarity  discriminant.  In  the 
circle  on  the  lower  left  are  the  earthquakes  (aftershocks  to  the  Mw  7.4  Suusamyr  thrust  earth¬ 
quake)  used  in  calibrating  the  discriminant  in  this  area.  The  middle  circle  encloses  the  earth¬ 
quakes  that  were  used  to  test  the  discriminant. 


tions  (calculated  by  forming  a  weighted  average  of  the  single  station  solutions  with  weighting  to  give 
preference  to  sonograms  with  higher  SNR)  appear  to  be  slightly  more  robust  (although  generally  only 
2  stations  voted  so  the  difference  is  slight). 

3.2  KNET: 

In  this  dataset  (Figures  16  and  18)  there  are  enough  events  to  permit  calibration  of  the  discriminant  and 
testing  on  an  independent  dataset.  The  calibration  of  the  discriminant,  using  quarry  blasts  north  of 
CHM  and  earthquakes  to  the  south  of  AML  (aftershocks  to  the  Mw  7.4  Suusamyr  thrust  earthquake), 
shows  that  single-station  solutions  are  highly  scattered  (Figure  19).  Calculating  network  solutions  we 
find  that  the  earthquake  population  is  tightly  grouped  and  little  crossover  with  the  quarry  blasts  exists 
(Figure  20).  Shown  in  this  figure  are  the  mean  and  +/-  1,2  and  3  standard  deviations. 
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KNET  events  —  single  station  solutions 


event  number 


source— receiver  ranges 


0  5  10  15  20  25  30  35  40  45 


event  number 

Figure  19:  Single  station  solutions  and  ranges.  Events  1  to  25  are  believed  to  be  quarry  blasts. 
The  remaining  events  are  aftershocks  to  the  Suusamyr  thrust  event. 

Using  earthquakes  located  near  station  AAK,  and  different  quarry  blasts  from  the  same  mine,  we  find 
that  the  network  solutions  are  essentially  the  same  (Figure  21).  The  standard  deviation  lines  are  from 
the  calibration  set  of  earthquakes  (Figure  20).  When  a  network  is  available  the  similarity  discriminant 
appears  to  be  quite  effective  and  robust  at  local  and  near-regional  ranges.  This  qualitative  assessment 
will  be  followed  shortly  by  a  more  rigorous  description  of  population  separation  and  misclassification 
probability  based  in  statistics. 
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full  network  solutions 
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Figure  20:  Network  solutions  using  the  same  (calibration)  events  considered  in  the  preceding 
figure.  On  the  bottom  we  display  the  number  of  stations  used  in  each  solution.  The  earthquake 
solutions  follow  a  normal  distribution.  The  horizontal  lines  represent  the  mean  and  +/-  1,2  and  3 
standard  deviations. 
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Figure  21;  Network  solutions  using  the  test  events  (see  Figure  20).  The  horizontal  lines  in  the 
top  figure  are  reproduced  from  Figure  22  (the  calibration  events). 


3.3  Outliers 

An  important  aspect  of  our  research  program  is  exploring  the  reasons  for  outliers.  In  section  3  the 
similarity  discriminant  produced  several  outliers.  In  this  section  we  consider  three  of  them.  Some 
outliers  can  be  explained  as  the  result  of  poor  station  coverage  (e.g.  BOSIET  event  70  which  was  re¬ 
corded  by  a  single  station).  Other  outliers  were  recorded  by  more  than  1  station  (e.g.  NRDC  events  1 
&  10  -  2  stations  and  KNET  event  26-9  stations).  Unless  our  a  priori  event  ID’s  are  incorrect  the 
similarity  algorithm  judged  the  NRDC  quarry  blast  to  be  an  explosion  and  the  KNET  earthquake  to  be 
a  quarry  blast. 

NRDC  events  1  &  10,  believed  to  be  quarry  blasts,  appear  to  have  generated  no  spectral  modulation  in 
the  band  from  0  to  50  Hz.  There  are  a  number  of  possible  explanations  -  the  delays  chosen  were  so 
short  that  no  modulations  were  observed  below  50  Hz.  In  theory  this  could  be  caused  by  delays  <  20 
ms  however  even  with  such  delays  we  might  expect  modulations  caused  by  source  finiteness.  Another 
alternative  is  that  due  to  shot  scatter  or  to  the  spatial  configuration  of  the  shots  there  was  not  enough 
shot  regularity  to  produce  modulations  at  any  frequency.  It  is  also  possible  that  these  events  were  not 
quarry  blasts  but  were  earthquakes,  due  to  natural  forces  or  related  to  the  mining  activity.  Unfortu¬ 
nately  without  “ground  truth”  information  we  can  only  speculate  as  to  the  cause  of  these  outliers. 

KNET  event  26  is  believed  to  be  an  earthquake  that  appears  to  have  a  significant  amount  of  time- 
independent  energy.  This  could  be  caused  by  near-source  resonance  however,  given  that  all  earth¬ 
quakes  in  the  group  used  to  calibrate  the  discriminant  are  tightly  grouped,  it  is  odd  that  the  other  events 
were  not  similarly  affected. 

4.  CONCLUSIONS  AND  FUTURE  PLANS 

4.1  Conclusions: 

It  appears  that  similarity  can  be  used  to  separate  quarry  blasts  from  instantaneous  explosions  and  earth¬ 
quakes.  When  used  to  discriminate  between  earthquakes  and  explosions  in  Kyrgyzstan  the  similarity 
discriminant  is  much  more  robust  when  a  network  is  available.  In  Kazakhstan  the  probability  of 
misclassifying  an  instantaneous  explosion  as  a  mining  event  is  higher  when  surface  stations  are  used  - 
probably  because  of  near-surface  resonance  at  some  sites. 

Outliers,  present  in  both  the  KNET  and  NRDC  analyses,  cannot  be  properly  explained  without  ground 
truth  information.  We  need  more  controlled  experiments  involving  quarry  blasts  and  calibration  explo¬ 
sions  with  ground  truth  information  for  all  events. 

4.2  Future  Plans: 

Although  this  result  comes  at  the  end  of  the  contract  period  it  is  preliminary  to  the  larger  effort  of 
developing  effective  means  to  discriminate  ripple-fired  events  from  instantaneous  explosions  and  earth¬ 
quakes.  This  will  be  continued  under  our  current  contract  F19628-95-K-0012.  We  are  currently  creat¬ 
ing  a  technique  that  merges  the  time-independence  and  component  independence  into  a  single  dis¬ 
criminant. 

In  addition  we  need  to  extend  this  research  by  analyzing  array  and  network  recordings  made  in  other 
regions  at  greater  ranges  to  test  robustness  and  transportability  and  determine  the  effective  range  of  the 
discriminants.  Wherever  possible  we  will  analyze  ground  truthed  data  to  better  understand  outliers. 
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should  be  produced  by  mines  that  use  short  (e.g.  20  ms)  delays  especially  when  delay  times  are  irregu¬ 
lar.  We  will  test  the  theory  that  source  finiteness  can  yield  a  low  frequency  time-independent  signature. 
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Appendix  A.  Sono 


1.  OVERVIEW 

The  purpose  of  the  sonogram  software  is  to  expand  time  series  into  2-dimensional  displays  which  show 
the  dependence  of  spectral  energy  on  time  and  frequency.  Our  research  into  the  time  -frequency  char¬ 
acteristics  of  single  and  multiple  chemical  explosions  detonated  and  recorded  in  Kazakhstan,  FSU 
(Hedlin  etal,  1989)  indicates  that  a  time  -independent  spectral  modulation  can  be  caused  by  ripple  - 
firing  (a  blasting  technique  in  which  a  number  of  closely  grouped  explosions  are  detonated  in  rapid 
sequence). . 

The  sonogram  software  expands  a  time  series  into  a  time  -frequency  display  by  using  a  multi  -taper 
spectral  estimation  procedure.  Successive  spectral  estimates  are  obtained  from  data  isolated  by  a  slid¬ 
ing  window.  Multi  -taper  spectral  estimation  has  been  used  because  of  the  proven  ability  of  this  ap¬ 
proach  to  yield  minimally  biased  spectral  estimates  from  short  time  series.  See  Park  et  al.  (1987)  for 
an  explanation  of  this  technique.  . 

Given  the  complex  nature  of  most  spectral  estimates  we  have  introduced  the  option  to  reduce  the 
sonograms  to  binary  form.  In  this  processing  each  spectral  estimate  is  converted  into  two  smoothed 
versions  which  are  sensitive  to  small  -  and  large  -scale  spectral  features.  These  versions  are  obtained 
by  convolving  the  spectral  estimate  with  a  short  and  long  boxcars  respectively.  . 

By  differencing  the  two  estimates  and  replacing  the  difference  at  each  frequency  with  a  +1  to  indicate 
locally  high  power  and  a  -1  to  indicate  the  opposite  we  convert  the  original  spectral  estimate  into  a 
binary  spectral  estimate.  Applying  this,  highly  processing  to  each  individual  spectral  estimate  in  the 
sonogram  yields  a  two  -dimensional  binary  pattern.  In  Hedlin  et  al.  (1989,  1990)  we  found  that  this 
processing  aided  discrimination  by  suppressing  minor  spectral  details  and  thus  enhancing  any  time  - 
independent  spectral  signature  that  might  be  present. 

Given  the  density  of  seismic  stations  in  some  areas  (due  to  regional  networks  or  arrays)  many  seismic 
events  are  redundantly  recorded.  As  a  result  it  is  possible  to  enhance  the  time  -frequency  spectral 
estimates  by  stacking.  As  discussed  in  Hedlin  et  al.  (1990),  stacking  can  be  performed  in  many  ways 
(e.g.  before  or  after  conversion  to  binary  form).  It  has  not  been  proven  which  approach  is  best.  As  a 
result  we  have  not  included  this  capability  in  this  version  of  the  sonogram  software.  Stacking  sonogram 
matrices  is  easily  performed  within  matlab. 

2.  APPLICATION 

In  the  figures  following  the  man  page  we  show  an  application  of  sono  to  events  recorded  by  the  NRDC 
network  in  Kazakhstan  (Figure  1).  In  Figure  2  we  show  the  time-series  and  sonogram  of  a  calibration 
explosion  (CH2  recorded  at  high  gain  by  the  vertical  component  surface  sensor  at  Bayanaul).  This 
sonogram  displays  no  organization  of  energy  into  time-independent  bands  (as  expected)  since  the  source 
is  an  instantaneous  explosion.  The  sonogram  displayed  in  Figure  3  shows  significant  banding.  This  is 
believed  to  be  the  result  of  ripple  firing  at  the  source.  In  Figure  4  to  6  we  schematically  show  the 
method  by  which  we  convert  spectral  sonograms  to  binary  form.  The  intent  of  this  conversion  is  to 
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suppress  large  scale  features  in  the  sonogram  and  enhance  small-scale  features  (like  time-indepen¬ 
dence)  that  might  indicate  ripple-firing. 
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NAME 

sono  -  calculates,  indexes,  and  stores  a  sonogram  (amplitudes  in  the  time-frequency  plane)  for  a 
seismogram. 

SYNOPSIS 

sono  -  returns  the  proper  "usage  line", 
sono  p2ir=parfile 

No  spaces  are  allowed  around  "="  on  the  command  file. 

ENVIRONMENT 

Sono  is  designed  to  operate  in  the  CSS  3.0  database  environment.  The  program  was  tested  with  both 
the  ORACLE  and  SYBASE  database  management  systems. 

For  users  of  the  SYBASE  DBMS,  the  following  environment  variables  should  be  set:  SYBASE, 
DSQUERY,  GDIHOME.  SYBASE  and  DSQUERY  specify  the  interface  with  the  DBMS,  and 
GDIHOME  is  set  to  the  top  directory  where  the  SAIC  Generic  Database  Interface  (GDI)  is  installed. 

For  example, 

setenv  SYBASE  /data/sybase/sybase492 
setenv  DSQUERY  SYBASE492 
setenv  PATH  $SYBASE/bin:$PATH 
setenv  GDIHOME  /data/dbjibs/iibgdi 

For  users  of  the  ORACLE  DBMS  the  following  environment  variables  should  be  set  ORACLE_HOME, 
ORACLE_SID. 

For  example, 

setenv  ORACLE_HOME  /data/oracle/6036 

setenv  ORACLE_SID  dev6036 

setenv  PATH  $ORACLE_HOME/bin:$PATH 

Users  must  have  a  valid  account,  password,  and  a  default  database  to  access  and  use.  This  database 
must  consist  at  least  of  the  wfdisc,  lastid  and  tfdisc  tables.  See  the  section  TFDISC  TABLE  below  for 
a  description  of  tfdisc. 


INPUT 

The  input  to  sono  follows  the  style  developed  at  SAIC  for  the  various  waveform  management  and  data¬ 
base  access  libraries.  User  input  is  specified  through  a  "parameter  file"  paifile  which  contains  an  ASCII 
description  of  user— specifiable  parameters.  A  parfile  consists  of  a  list  of  parameters  and  values  of  the 
form  parameter=value,  as  required  by  the  program.  No  spaces  are  allowed  around  "="  in  the  file.  All 
parameter  names  should  be  spelled  exactly  as  they  appear  below  (see  section  PARAMETERS).  All 
parameter  values  that  are  character  strings  are  case-sensitive  unless  noted  otherwise. 

Parameters  can  appear  in  any  order.  If  any  parameter  needed  for  proper  execution  is  missing,  then  the 
program  exits  with  an  error  message  informing  the  user  of  the  missing  information.  Comments  are 
allowed  anywhere  in  the  file,  but  not  on  the  same  line  as  a  parameter.  Comment  should  not  start  with  a 
valid  parameter  name.  For  stylistic  reasons,  it  is  preferable  to  adopt  a  single  convention  for  the  com¬ 
ments  syntax. 

The  parameters  in  the  file  can  be  divided  into  two  sections  according  to  their  use  by  the  program.  The 
first  section  provides  parameters  for  the  database  connection  and  access.  The  second  section  includes 
parameters  which  control  the  sonogram  algorithm  and  calculation. 
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Example  of  parfile  : 

#  The  following  line  is  just  a  comment 
###############  for  DB  connection  ############## 
in_db=nrdc/nrdc 

account=''nrdc" 

password="nrdc'’ 

vendor=”sybase” 

#  this  is  a  comment  line 
staJist="BAY” 
chanJist=”HSSE" 

s  tart_time=54 8073338.0 

#  if  preferred,  use  end_time  instead  of 

#  duration,  but  not  both 
duration=35.0 
extension_time=360 
join=2 
wfd_to_wfm 

#  The  following  line  is  just  a  comment 
############  for  sonogram  calculation  ############# 
BinSono=true 

NarrowBoxCar=2.5 

BroadBoxCar=5,0 

OutputFormat=matIab 

OutputDir=/data/sonos 

OutputTable=t_tfdisc 

#  The  next  two  parameters  can  have  values 

#  displacement,  velocity,  or  acceleration 
DataTypeIn=velocity 
DataTypeOut=velocity 
WindowSlide=20.0 
WTimeLength=2.0 
LowFreq=10.0 

HighFreq=35.0 

TimeBandwidth=4.0 

################################################### 


PARAMETERS 

FIRST  SECTION 

All  parameters  in  this  section  are  not  used  by  the  sono  program  itself,  but  are  passed  to  the  Libgdi  and 
Libwfm  libraries.  Sono  only  provides  minimal  checking  of  these  parameters,  and  relies  on  the  error 
handling  capabilities  of  the  libraries  to  deal  with  improper  input. 

account=char 

User  name  as  a  login  account  in  the  database. 

Use  double  quotes  around  the  character  string. 

Example  for  SYBASE  :  account="nrdc" 

Example  for  ORACLE  :  account=’'nrdc/nrdc@t:nrd:dev71 3" 

chan  list=char 

Contains  the  name  of  the  channel  from  which  a  sonogram  is  to  be  calculated. 
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Use  double  quotes  around  the  character  string. 

Example:  chan_list="HSSE'’ 

See  also  libwfm(3) 

duration=double 

Desired  length  (in  seconds)  of  the  seismogram  from  which  to  calculate  the  sonogram,  start¬ 
ing  from  the  time  specified  as  start_time. 

Either  duration  or  end_tinie  should  be  specified, 
but  not  both. 

Example:  duration=35.0 
See  also  libwfm(3) 

end_time=double 

Desired  end_time  for  the  seismogram  and  sonogram  (in  seconds). 

Either  duration  or  end_tinie  should  be  specified, 
but  not  both. 

Example:  end_time=548073352.0 
See  also  libwfm(3) 

extension_time=double 

Number  of  seconds  to  be  subtracted  from  start_time  and  added  to  end_time  before  search¬ 
ing  the  wfdisc  table  (in  order  to  make  the  search  more  efficient.)  A  good  default  value  to 
use  should  be  of  the  order  of  minutes. 

Example:  extension_time=360 
See  also  Iibwfm(3) 

in_db=char/char 

Input  database  account  and  password 

Example  for  SYBASE:  in_db=nrdc/nrdc 

Example  for  ORACLE:  in_db=nrdc/nrdc@t:nrd:dev713 

See  also  libwfm(3) 


join=:int 

Can  be  set  to  0,  1,  2.  The  default  value  of  0  does  not  allow  adjacent  waveforms  to  be  con¬ 
catenated.  A  value  of  1  means  that  input  waveform  data  adjacent  in  time  are  joined  and 
indexed  through  one  output  wfdisc.  Time  gaps  will  be  filled  with  values  defined  by  the 
default  fill  value  used  in  the  Libwfm  library,  namely  -1,0.  A  value  of  2  requests  that  a 
conditional  join  be  performed,  i.e.  adjacent  waveform  data  will  be  joined  provided  the  time 
gap  between  them  is  less  than  gapepsilon,  which  is  defaulted  to  the  inverse  of  the  sampling 
rate  as  per  the  Libwfm  convention.  A  value  of  2  is  preferred  unless  the  user  is  sure  that 
another  choice  is  better. 

Example:  join=2 
See  also  Iibwfm(3) 

password=char 

User  password  for  the  login  account  in  the  database. 

Use  double  quotes  around  the  character  string. 

Example:  password='’nrdc” 

stalist=char 

Contains  the  name  of  the  station  of  interest. 

Use  double  quotes  around  the  character  string. 

Example:  sta_list="BAY" 
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See  also  Iibwfm(3) 
start_tinie=double 

The  beginning  time  of  interest  for  the  calculation  of  the  sonogram  (in  seconds). 
This  is  an  epoch  time. 

Example:  start_time=548073338.0 

vendor=char 

Vendor  name  for  the  database. 

Tested  databases:  oracle,  Sybase. 

Use  double  quotes  around  the  character  string. 

Example:  vendor=''sybase” 


verbose 

Flag.  When  used,  the  information  about  the  libwfm  run  is  printed  to  stdout.  To  suppress  the 
printing  don't  set  this  flag. 

Example:  verbose 
See  also  Iibwfin(3) 

wfdtowfm 

Flag.  Indicates  that  the  created  wfmem  information  must  be  stored  in  core  memory. 
Required. 

Example:  wfd_to_wfm 
See  also  libwf)m(3) 

SECOND  SECTION 

For  parameters  in  this  section  the  sono  program  provides  an  error  check  for  their  names  and  values,  and 
issues  a  diagnostic  error  message. 

B  inSono=boolean 

Boolean  flag,  takes  values  true  or  false  case  insensitive.  When  set  to  true  the  sonogram  is 
converted  to  a  matrix  of  values  of  +1  or  -1,  (labeled  binary  sonogram  by  Hedlin  et  al., 
1989)  according  an  the  algorithm  outlined  in  the  DESCRIPTION  section  below.  Otherwise 
sonogram  values  are  left  as  double  precision  numbers. 

Example:  BinSono=false 

NarrowBoxCar=double 

This  parameter  is  required  only  if  BinSono=true.  Indicates  a  length  (in  seconds)  of  a  nar¬ 
row  boxcar  function  used  to  smooth  each  spectrum  prior  to  conversion  to  a  binary  sono¬ 
gram 

Example:  NarrowBoxCar=2.5 

BroadBoxCar=double 

This  parameter  is  required  only  if  BinSono=true.  Indicates  a  length  (in  seconds)  of  a 
broad  boxcar  function  used  to  smooth  each  spectrum  prior  to  conversion  to  a  binary  sono¬ 
gram 

Example:  BroadBoxCar=5.0 

OutputFormat=char 

Specifies  the  output  format  of  the  sonogram.  Case  insensitive.  Available  formats  are 
matlab,  binary,  ascii.  See  section  OUTPUT  of  this  man  page  for  formats  details. 

Example:  OutputFormat=matlab 
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OutputDir=char 

Full  path  name  of  a  directory  where  the  resulting  sonogram  is  placed. 

Example:  OutputDir=/data/sonograms/l 995/binary _sonos 

Ou  tputTable=char 

Name  of  the  database  table  which  is  updated  with  information  about  the  calculated  sono¬ 
gram.  For  this  purpose  a  new  table  tfdisc  (according  to  CSS  3.0  convention)  was  designed. 
See  section  TFDISC  TABLE  of  this  man  page. 

Example:  OutputTable=tfdisc 

DataTypeIn=char 

Flag  indicating  what  type  of  data  is  indexed  by  the  wfdisc  entry.  This  flag  is  used  to  indi¬ 
cate  the  natural  units  of  the  detector. 

Allowed  values  are  displacement,  velocity  and  acceleration. 

Example:  DataTypeIn=velocity 

DataTypeOut=char 

Flag  indicating  what  type  of  data  will  be  written  to  the  resulting  sonogram.  Allowed  values 
are  displacement,  velocity  and  acceleration. 

Example:  DataTypeOut=velocity 

WindowSIide=double 

Indicates  by  how  much  to  shift  the  time  window  for  successive  spectra  in  the  sonogram,  as 
a  percentage  of  window  length. 

Example:  WindowSlide=20.0 

WTimeLength=double 

Size  of  the  time  window  (in  seconds)  used  to  compute  the  successive  spectra  in  the  sono¬ 
gram. 

Example:  WTimeLength=2.0 

LowF  req=doubIe 

Lower  limit  of  the  frequency  range  of  interest.  Negative  input  is  reset  to  0.0 
Example:  LowFreq=10.0 

HighFreq=double 

Upper  limit  of  the  frequency  range  of  interest. 

Example:  HighFreq=35.0 

TimeBandwidth=4.0 

Time  bandwidth  used  ... 

Example:  TimeBandwidth=4.0 


OUTPUT 

The  main  output  of  sono  is  a  sonogram,  in  one  of  several  user  specifiable  formats.  In  addition,  the  pro¬ 
gram  updates  the  database  table  tfdisc  with  appropriate  indexing  (or  header)  information.  (See  section 
TFDISC  TABLE  below  for  details  on  the  tfdisc  contents.) 

The  file  name  is  constructed  from  the  station  and  channel  names,  and  waveform  date  and  start  time 
used  in  the  sonogram  calculation.  Our  present  convention  produces  file  names  of  the  form  . 
sta  chan  yyyydddhhmmss.tf,  where 
sta  -  station  name 
chan  -  channel  name 
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yyyy  “  y^^^ 

ddd  -  Julian  day 
hh  -  hours 
mm  ~  minutes 
ss  -  seconds 

Currently  one  of  three  different  formats  may  be  selected:  matlab,  binary,  ascii.  The  program  architec¬ 
ture  allows  easy  expansion  for  additional  user’s  formats.  If  such  need  arises  the  appropriate  additions 
should  be  made  to  the  functions  map_format_type(),  unmap_format_type()  in  the  file  mapping.c,  and  to 
the  function  write_sono()  in  the  file  write_sono.c. 

In  the  current  formats  Nrows  is  the  number  of  frequencies  and  Ncolumns  is  the  number  of  time  win¬ 
dows  used  for  the  calculations.  Each  data  point  in  the  matrix  is  an  8-byte  double  precision  number  (in 
case  BinSono=false)  or  4-byte  signed  integer  (in  case  BinSono=true) 

matlab  - 

The  sonogram  file  structure  is  that  of  a  Level  1.0  MAT-file.  The  file  consists  of  a  2-part 
header  followed  by  the  matrix  of  sonogram  data  points. 

The  first  part  of  the  header  is  a  standard  20-byte  matlab  header  which  consists  of  5  four-byte 
long  integers: 

1.  type  - 

type  flag.  An  integer  whose  4  decimal  digits  encode  storage  information.  The  first 
digit  indicates  the  numeric  format  of  binary  numbers  on  the  computer  that  wrote  the 
file.  Sono  uses  1,  indicating  IEEE  Big  Endian  (Macintosh,  SPARC,  Apollo,  SGI, 
HP9000/300).  Other  values  can  be  obtained  from  the  matlab  documentation  if  the 
software  is  to  be  ported  to  a  computer  with  a  different  architecture.  The  second  digit 
is  always  0.  The  third  digit  can  take  one  of  two  values:  0  -  which  means  that  data  are 
stored  as  double  precision,  this  is  when  no  conversion  to  a  binary  sonogram  was 
specified.  2  -  data  are  stored  as  signed  integers,  this  is  when  the  conversion  to  a 
binary  sonogram  was  requested.  The  fourth  digit  is  0  which  merely  indicates  that  the 
matrix  is  numeric  instead  of  consisting  of  text  items.  To  summarize,  the  sono  pro¬ 
gram  uses  two  types,  both  defined  in  sono.h. 
type=1020  -  when  writing  double  precision  values  and 
type:=1000  -  when  writing  signed  integer  values. 

2.  Nrows  - 

number  of  rows  in  the  matrix. 

3.  Ncolumns 

number  of  columns  in  the  matrix. 

4.  Iflag  “ 

imaginary  flag.  Always  0  since  the  resulting  sonogram  has  no  imaginary  part,  only 
real. 

5.  NameLength  - 

number  of  bytes  to  hold  the  matrix  name  plus  1  byte  for  the  NULL  character. 

The  second  part  of  the  header  contains  a  matrix  name  +  1  for  the  NULL  character.  The  matrix 
name  is  defined  in  sono.h  file  by  MATNAME.  Total  size  of  the  header  is  <20  +  namelength  + 
1>  bytes.  This  name  is  of  no  importance  to  the  analyst,  but  may  be  used  and/or  modified  by 
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MATLAB  users  to  manipulate  the  matrix. 

Sonogram  data  points  are  stored  in  the  file  immediately  after  the  header  as  a  matrix.  The 
matrix  consists  of  Nrows  *  Ncolumns  numbers  which  are  stored  column-wise  in  binary  format. 


binary  - 

The  output  file  consists  of  the  sonogram  points  presented  as  a  matrix,  where  Nrows  « 
Ncolumns  numbers  are  stored  column-wise.  Each  number  requires  8  bytes  for  double  precision 
or  4  bytes  for  signed  integer.  With  this  choice  the  output  file  contains  no  header  and  it  is 
assumed  that  the  proper  indexing  is  done  exclusively  through  the  DBMS  entry  in  the  tfdisc 
table. 

ascii  -  The  output  file  consists  of  the  sonogram  points  presented  as  a  matrix,  where  Nrows 

Ncolumns  numbers  are  stored  column-wise  and  written  in  the  file  one  per  line  in  the  format 
"%15.7E\n"  for  double  precision  or  "%d"  for  signed  integer.  With  this  choice  the  output  file 
contains  no  header  and  it  is  assumed  that  the  preper  indexing  is  done  exclusively  through  the 
DBMS  entry  in  the  tfdisc  table. 

DESCRIPTION 

The  purpose  of  the  sonogram  software  is  to  expand  time  series  into  2-dimensional  displays  which  show 
the  dependence  of  spectral  energy  on  time  and  frequency.  Our  research  into  the  time— frequency  charac¬ 
teristics  of  single  and  multiple  chemical  explosions  detonated  and  recorded  in  Kazakhstan,  FSU  (Hedlin 
et  ah,  1989)  indicates  that  a  time-independent  spectral  modulation  can  be  caused  by  ripple-firing  (a 
blasting  technique  in  which  a  number  of  closely  grouped  explosions  are  detonated  in  rapid  sequence). 

The  sonogram  software  expands  a  time  series  into  a  time-frequency  display  by  using  a  multi-taper 
spectral  estimation  procedure.  Successive  spectral  estimates  are  obtained  from  data  isolated  by  a  sliding 
window.  Multi-taper  spectral  estimation  has  been  used  because  of  the  proven  ability  of  this  approach 
to  yield  minimally  biased  spectral  estimates  from  short  time  series.  See  Park  et  al.  (1987)  for  an  expla¬ 
nation  of  this  technique. 


Given  the  complex  nature  of  most  spectral  estimates  we  have  introduced  the  option  to  reduce  the  sono¬ 
grams  to  binary  form.  In  this  processing  each  spectral  estimate  is  converted  into  two  smoothed  ver¬ 
sions  which  are  sensitive  to  small—  and  large-scale  spectral  features.  These  versions  are  obtained  by 
convolving  the  spectral  estimate  with  a  short  and  long  boxcars  respectively. 

By  differencing  the  two  estimates  and  replacing  the  difference  at  each  frequency  with  a  +1  to  indicate 
locally  high  power  and  a  —1  to  indicate  the  opposite  we  convert  the  original  spectral  estimate  into  a 
binary  spectral  estimate.  Applying  this,  highly  non-linear,  processing  to  each  individual  spectral  esti¬ 
mate  in  the  sonogram  yields  a  two-dimensional  binary  pattern.  In  Hedlin  et  al.  (1989,  1990)  we  found 
that  this  processing  aided  discrimination  by  suppressing  minor  spectral  details  and  thus  enhancing  any 
time-independent  spectral  signature  that  might  be  present. 

Given  the  density  of  seismic  stations  in  some  areas  (due  to  regional  networks  or  arrays)  many  seismic 
events  are  redundantly  recorded.  As  a  result  it  is  possible  to  enhance  the  time- frequency  spectral  esti¬ 
mates  by  stacking.  As  discussed  in  Hedlin  et  al.  (1990),  stacking  can  be  performed  in  many  ways  (e.g. 
before  or  after  conversion  to  binary  form).  It  has  not  been  proven  which  approach  is  best.  As  a  result 
we  have  not  included  this  capability  in  this  version  of  the  sonogram  software.  Stacking  sonogram 
matrices  is  easily  performed  within  matlab. 


FILES 

sono.h  include  file  for  C— source  files 

sono  f77.h  include  file  for  FORTRAN-sorces  files 
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db  tfdisc.h  definition  of  a  tfdisc  structure 

DIAGNOSTICS 

The  routines  produce  error  messages  on  stderr.  Most  messages  report  the  function  that  was  active  when 
the  error  have  occurred. 

LIMITATIONS 

This  code  does  not  provide  for  interactive  selection  of  time  windows,  etc.  This  choice  is  deliberate, 
because  a  graphical  interface  should  be  provided  by  existing  modules  at  the  time  of  integration  with 
other  software. 

Pre-processing  steps  such  as  frequency-filtering  of  the  waveform  data  prior  to  sonogram  calculation  are 
assumed  to  have  been  accomplished  using  other  modules.  If  necessary  the  code  could  be  altered  to 
accommodate  such  tasks.  This  would  be  best  done  at  integration  time,  since  the  appropriate  modules 
would  not  have  to  be  duplicated. 

Similarly,  stacking  of  the  resulting  sonograms,  or  beamforming  prior  to  sonogram  calculation  is  not  pro¬ 
vided  as  part  of  this  code.  Such  steps  are  application  and  user-dependent,  and  are  best  accomplished 
within  other  modules,  or  even  in  the  context  of  other  general-purpose  packages  such  as  MATLAB. 

BUGS 

No  known  bugs  at  this  point. 

Report  bugs  to  nwilliams@ucsd.edu 

LANGUAGE  INTERFACE 
C,  FORTRAN 

SEE  ALSO 

Hedlin,  M.  A.  H.,  Minster,  J.  B.,  &  Orcutt,  J.  A.,  1989.  The  time-frequency  characteristics  of  quarry 
blasts  and  calibration  explosions  recorded  in  Kazakhstan,  USSR,  Geophysical  Journal  International,  99, 
109-121. 

Hedlin,  M.A.H.,  Minster,  J.B.  &  Orcutt,  J.A.,  1990.  An  automatic  means  to  discriminate  between  earth¬ 
quakes  and  quarry  blasts,  BulL  SeismoL  Soc.  Am.,  80,  2143-2160. 

Park,  J.,  Lindberg,  C.R.  and  Vernon,  F.L.,  1987.  Multitaper  spectral  analysis  of  high-frequency  seismo¬ 
grams,  J.  Geophys.  Res.,  92,  12675-12684. 

Iibwfm(3),  libwfdisc(3) 

AUTHOR 

Michael  A.H.  Hedlin,  UCSD/SIO/IGPP,  San  Diego,  CA  92093-0225 
Nadya  P.  Williams,  UCSD/SIO/IGPP,  San  Diego,  CA  92093-0225 


APPENDIX 


TFDISC  TABLE 

The  following  contains  the  description  of  a  new  table  intended  as  an  extension  to  the  CSS  3.0  database 
schema,  to  describe  sonograms,  sono  updates  this  relation  automatically  when  a  sonogram  is  calcu¬ 
lated.  Sono  supports  either  the  ORACLE  or  SYBASE  DBMS,  through  the  SAIC  Generic  Database 
Interface  (GDI).  The  new  table  is  meant  to  describe  a  general  time-frequency  representation  of  a  time 
series. 


Relation  :  tfdisc 

Description  :  sonogram  file  header  and  descriptive  information 


Columns  Description: 

1  -  Attribute  NAme 
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2  -  Field  Number 

3  ~  Storage  Type 

4  -  External  Format 

5  “  Character  Position 

6  ~  Attribute  Description 


1 

2 

3 

4 

5 

6 

sta 

1 

c6 

a6 

1-6 

station 

chan 

2 

c8 

aS 

8-15 

channel 

time 

3 

f8 

fl7.5 

17-33 

epoch  time 

tfid 

4 

i4 

i8 

35^2 

sonogram  id 

jdate 

5 

i4 

i8 

44-51 

Julian  date 

dsono 

6 

f8 

f9.2 

53-61 

duration  of  sonogram 

dwind 

7 

f8 

f6.2 

63-68 

duration  of  each  window 

pslide 

8 

f8 

f5.1 

70-74 

offset  between  windows 

nt 

9 

i4 

i4 

76-79 

number  of  time  windows 

fmin 

10 

f8 

f7.4 

81-87 

minimum  frequency 

fmax 

11 

f8 

f7.4 

89-95 

maximum  frequency 

nf 

12 

i4 

i4 

97-100 

number  of  frequencies 

tbprod 

13 

f8 

f5.2 

102-106 

time-bandwidth  product 

neigspect 

14 

i4 

i4 

108-111 

number  of  eigenspectra 

spectype 

15 

c2 

a2 

113-114 

spectrum  type 

bbroad 

16 

f8 

f4.1 

116-119 

broad  boxcar  filter 

bn  arrow 

17 

f8 

f4.1 

121-124 

narrow  boxcar  filter 

foff 

18 

i4 

i4 

126-129 

byte  offset 

dir 

19 

c64 

a64 

131-194 

directory 

dfile 

20 

c32 

a32 

196-227 

tf  file 

tfformat 

21 

c9 

a9 

229-237 

tf  file  format 

progver 

22 

cl5 

al5 

239-253 

program  identifier 

Iddate 

23 

date 

al7 

255-271 

load  date 

TFDISC  TABLE  ATTRIBUTES 


Name: 

Description: 

ORACLE: 

SYBASE: 

NA  Value: 

Range: 

sta 

Station  code.  This  is  a  common  code-name  of 
seismic  observatory.  Usually  3  or  4  characters 
are  used. 

VARCHAR(6) 

VARCHAR(6) 

A  valid  entry  is  required. 

Any  upper  case  string  up  to  6  characters  long. 

Name: 

chan 

Description: 

Channel  identifier.  This  is  an  eight-character 
code  which,  taken  together  with  sta  and  time, 
uniquely  identifies  the  source  of  the  seismic  data. 

ORACLE: 

VARCHAR(8) 

SYBASE: 

VARCHAR(8) 

NA  Value: 

A  valid  entry  is  required. 

Range: 

Any  upper  case  string  up  to  6  characters  long. 
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Name: 

Description: 

ORACLE: 

SYBASE: 

NA  Value: 

Units: 

time 

Epoch  time.  Given  in  seconds  and  fractions  of  a 
second  since  hour  0  January  1,  1970.  Stored  in  a 
double  precision  floating  point  number. 

FLOAT(53) 

FLOAT 

A  valid  entry  is  required. 

Seconds 

Name: 

tfid 

Description: 

Sonogram  identifier.  The  key  field  is  a  unique 
identifier  for  a  sonogram  calculated  for  a  certain 
segment  of  waveform  data. 

ORACLE: 

NUMBER(8) 

SYBASE: 

INT 

NA  Value: 

A  valid  entry  is  required. 

Range: 

tfid  >  0 

Name: 

jdate 

Description: 

Julian  date.  This  attribute  is  the  date  of  an 

onset  of  a  P-wave  for  a  certain  segment  of  waveform 

data  for  which  the  sonogram  was  calculated. 

ORACLE: 

NUMBER(8) 

SYBASE: 

INT 

NA  Value: 

A  valid  entry  is  required. 

Range: 

Julian  dates  of  the  form  yyyyddd.  Must  be  consistent 
with  the  accompanying  time  attribute. 

Name: 

dsono 

Description: 

Duration  of  the  waveform  to  use  for  the 
sonogram  calculation. 

ORACLE: 

FLOAT(24) 

SYBASE: 

FLOAT 

NA  Value: 

A  valid  entry  is  required. 

Range: 

A  positive  integer.  Generally  a  value  less  than 

35  may  not  give  a  good  sonogram. 

Units: 

Seconds 

Name: 

dwind 

Description: 

Duration  of  a  time  window  during  the  sonogram 
calculation. 

ORACLE: 

FLOAT(24) 

SYBASE: 

FLOAT 

NA  Value: 

A  valid  entry  is  required. 

Range: 

A  positive  value  of  the  order  of  a  few  seconds. 

Units: 

Seconds 

Name: 

pslide 

Description: 

Offset  between  windows  as  a  percent  of  a 
window  duration. 

ORACLE: 

FLOAT(24) 

SYBASE: 

FLOAT 
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NA  Value; 

Range: 

Units: 

valid  value  is  required 

0  <  pslide  <  dwind 
percent 

Name: 

nt 

Description: 

Number  of  the  time  windows  used  for  the 
calculation  of  the  sonogram. 

ORACLE: 

FLOAT(24) 

SYBASE: 

FLOAT 

NA  Value: 

A  valid  entry  is  required. 

Range: 

nt  >  1 

Name: 

fmin 

Description: 

Minimum  frequency  used  for  the  calculation  of 
the  sonogram 

ORACLE: 

FLOAT(24) 

SYBASE: 

FLOAT 

NA  Value: 

A  valid  entry  is  required 

Range: 

0.0  <  fmin  <  fmax 

Units: 

Hertz 

Name: 

fmax 

Description: 

Maximum  frequency  used  for  the  calculation  of 
the  sonogram 

ORACLE: 

FLOAT(24) 

SYBASE: 

FLOAT 

NA  Value: 

A  valid  entry  is  required. 

Range: 

fmin  <  fmax  <  Nyquist  frequency 

Units: 

Hertz 

Name: 

nf 

Description: 

Number  of  frequencies  used  for  the  calculation 
of  the  sonogram 

ORACLE: 

FLOAT(24) 

SYBASE: 

FLOAT 

NA  Value: 

A  valid  entry  is  required 

Range: 

Is  calculated  by  the  program  from  dwind,  fmax, 
fmin  and,  if  exist,  bbroad  and  bnarrow. 

Name: 

tbprod 

Description: 

time-bandwidth  product 

ORACLE: 

FLOAT(24) 

SYBASE: 

FLOAT 

NA  Value: 

A  valid  value  is  required. 

Range: 

tbprod  >  0 

Name: 

neigspect 

Description: 

Number  of  Eigen  spectra. 

ORACLE: 

FLOAT(24) 

SYBASE: 

FLOAT 

NA  Value: 

A  valid  value  is  required 

Range: 

Is  calculated  by  the  program  from  tbprod 
as  <2*tbprod  -  l> 
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Name: 

Description: 


ORACLE: 
SYBASE: 
NA  Value: 
Range: 


spectype 

Flag  indicating  what  type  of  data  is  written 
to  the  resulting  sonogram,  and  whether  the  conversion 
to  the  binary  sonogram  was  performed 
VARCHAR(2) 

VARCHAR(2) 

A  valid  entry  is  required 

{  bd  I  bv  I  ba  I  d  I  V  I  a  },  lower  case 

These  flags  stand  for  : 

bd  ”  data  type  displacement, 

conversion  to  binary  sonogram 
bv  ~  data  type  velocity, 

conversion  to  binary  sonogram 
ba  -  data  type  acceleration, 

conversion  to  binary  sonogram 
d  -  data  type  displacement, 

no  conversion  to  binary  sonogram 
V  -  data  type  velocity, 

no  conversion  to  binary  sonogram 
a  ~  data  type  acceleration, 

no  conversion  to  binary  sonogram 


Name: 

Description: 

ORACLE: 

SYBASE: 

NA  Value: 

Range: 

Units: 

bbroad 

This  is  a  length  of  a  broad  boxcar  function 

used  to  smooth  each  spectrum  prior  to  conversion  to 

binary  sonogram 

FLOAT(24) 

FLOAT 

A  valid  entry  is  required 
bnarrow  <  bbroad  <  (fmax  -  fmin) 

If  no  BinSono  was  specified  then  this  value  is  0.0. 

Hertz 

Name: 

bnarrow 

Description: 

This  is  a  length  of  a  narrow  boxcar  function 
used  to  smooth  each  spectrum  prior  to  conversion 

to  binary  sonogram 

ORACLE: 

SYBASE:FLOAT 

FLOAT(24) 

NA  Value: 

A  valid  entry  is  required 

Range: 

0.0  <  bnarrow  <  bbroad 

If  no  BinSono  was  specified  then  this  value  is  0.0. 

Units: 

Hertz 

Name: 

foff 

Description: 

File  offset.  This  is  the  byte  offset 

(header  size)  of  the  data  segment  within  a  sonogram. 

ORACLE: 

NUMBER(8) 

SYBASE: 

INT 

NA  Value: 

A  valid  value  is  required 

Range: 

foff  >  0 
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Name: 

Description: 


ORACLE: 
SYBASE: 
NA  Value: 
Range: 


dir 

Directory.  A  directory  part  of  the  a  path  name 
to  the  location  of  the  created  sonogram.  A  relative 
path  or  a  may  be  used. 

VARCHAR(64) 

VARCHAR(64) 

A  valid  value  is  required. 

Any  string  up  to  64  characters  long,  where  each  character 
is  a  legal  character  for  a  directory  name  under  UNIX. 


Name: 

Description: 

ORACLE: 

SYBASE: 

NA  Value: 

Range: 

dfile 

File  name.  This  is  the  name  of  the  sonogram  file. 
VARCHAR(32) 

VARCHAR(32) 

A  valid  entry  is  required. 

Any  string  up  to  32  characters  long,  where  each  character 
is  a  legal  character  for  a  file  name  under  UNIX. 

Name: 

tfformat 

Description: 

Output  format  name.  A  name  of  a  format  used 
to  write  the  calculated  sonogram. 

ORACLE: 

VARCHAR(IO) 

SYBASE: 

VARCHAR(IO) 

NA  Value: 

A  valid  entry  is  required. 

Range: 

{  matlab  1  binary  1  ascii  },  lower  case 

Name: 

progver 

Description: 

Identifier  for  the  program’s  version. 

ORACLE: 

VARCHAR(IO) 

SYBASE: 

VARCHAR(IO) 

NA  Value: 

A  valid  entry  is  required 

Range: 

Any  string  of  10  characters  which  specifies 
the  version  and  the  name 

Name: 

Iddate 

Description: 

Load  date  This  is  the  date  and  time  when  the 
record  was  inserted  into  the  table. 

ORACLE: 

DATE 

SYBASE: 

DATETIME 

NA  Value: 

A  valid  entry  is  required 

Range: 

Any  valid  date 
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Figure  Al.  Geographic  locations  of  receivers  and  events  in  the  1989  NRDC  dataset 


Figure  A2.  Seismogram  resulting  from  the  single-chemical  explosion  Chemex  2  recorded 
at  a  range  of  157  km  by  the  vertical  component  surface  seismometer  at  Bayanaui 
and  corresponding  sonogram. 


Figure  A3.  Seismogram  resulting  from  event  c  recorded  at  a  range  of  264  km  by  the  vertical 
component  surface  seismometer  at  Bayanaul  and  corresponding  sonogram.  This  event  is 
believed  to  be  ripple  fired.  The  sonogram  matrix  contains  a  clear  time-independent  spectral 
signature. 


frequency  (Hz) 


Figure  A4.  Original  spectral  estimate  (solid)  and  two  filtered  versions,  one  relatively 
unsmoothed  (fine  dash)  and  the  other  heavily  smoothed  (coarse  dash).  This  figure  is 
intended  to  illustrate  the  means  by  which  we  reduce  each  spectral  estimate  to  a  binary 
spectral  estimate.  In  the  following  figures  regions  of  locally  high  power  are  represented  by 
a  +1  (white).  Regions  deficient  in  power  are  represented  by  a  -1  (black). 
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time(s) 
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Figure  A6.  Vertical  component  binary  sonogram  matrix  corresponding  to  event  c 
recorded  at  Bayanaul  (see  Figure  3).  Note  the  clear  time-independent  pattern. 
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